{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hilpisch.com/tpq_logo.png\" alt=\"The Python Quants\" width=\"35%\" align=\"right\" border=\"0\"><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Listed Volatility and Variance Derivatives\n",
    "\n",
    "**Wiley Finance (2017)**\n",
    "\n",
    "Dr. Yves J. Hilpisch | The Python Quants GmbH\n",
    "\n",
    "http://tpq.io | [@dyjh](http://twitter.com/dyjh) | http://books.tpq.io\n",
    "\n",
    "<img src=\"http://hilpisch.com/../images/lvvd_cover.png\" alt=\"Derivatives Analytics with Python\" width=\"30%\" align=\"left\" border=\"0\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Modeling of the VSTOXX Index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter is a bit different in style compared to the other chapters about the VSTOXX volatility index and related derivatives. It introduces another, parsimonious model, which is called square-root jump diffusion (SRJD) to model the VSTOXX volatility index. The model, which is essentially an extension of the Gruenbichler and Longstaff (1996) model as analyzed in the previous chapter, is capable of reproducing prices of European options written on the VSTOXX reasonably well.\n",
    "\n",
    "Two major enhancements characterize the SRJD model:\n",
    "\n",
    "* **term structure**: it allows to capture the term structure as observed in the prices of futures on the VSTOXX index\n",
    "* **jump component**: including a jump component allows to better replicate option prices in the short term\n",
    "\n",
    "Adding these two components makes a market-consistent calibration of the model to a comprehensive set of European options on the VSTOXX index possible. For similar analyses and modeling approaches including jumps for the volatility index refer, for example, to Psychoyios (2005), Sepp (2008) or Psychoyios et al. (2010)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Market Quotes for Call Options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we introduce the model, let us set the stage by revisiting market quotes for European call options written on the VSTOXX volatility index. First, we import the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "np.set_printoptions(suppress=True)\n",
    "path = './data/'\n",
    "h5 = pd.HDFStore(path + 'vstoxx_data_31032014.h5', 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have both options data and futures data stored in this file with quotes from 31. March 2014."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Option market quotes are what we are concerned with for the moment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "option_quotes = h5['options_data']\n",
    "option_quotes['DATE'] = pd.to_datetime(option_quotes['DATE'])\n",
    "option_quotes['MATURITY'] = pd.to_datetime(option_quotes['MATURITY'])\n",
    "option_quotes.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "option_quotes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At any given point in time, there are options on the VSTOXX available for eight maturities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "mats = sorted(set(option_quotes['MATURITY']))\n",
    "mats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spot level of the VSTOXX index on 31. March 2014 has been 17.6639. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v0 = 17.6639"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In what follows, we only want to plot call option market quotes which are not too far in-the-money nor out-of-the-money."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "tol = 0.4\n",
    "to_plot = option_quotes[(option_quotes['STRIKE'] > (1 - tol) * v0)\n",
    "                      & (option_quotes['STRIKE'] < (1 + tol) * v0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the VSTOXX European call option quotes which fulfill the requirements. The goal of this chapter is to replicate \"all these option quotes\" as good as possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from pylab import mpl, plt\n",
    "plt.style.use('seaborn')\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "markers = ['.', 'o', '^', 'v', 'x', 'D', 'd', '>', '<']\n",
    "plt.figure(figsize=(10, 6));\n",
    "for i, mat in enumerate(mats):\n",
    "    strikes = to_plot[(to_plot['MATURITY'] == mat)]['STRIKE']\n",
    "    prices = to_plot[(to_plot['MATURITY'] == mat)]['PRICE']\n",
    "    plt.plot(strikes, prices, 'b%s' % markers[i], label=str(mat)[:10])\n",
    "plt.legend();\n",
    "plt.xlabel('strike');\n",
    "plt.ylabel('option quote');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">VSTOXX European call option quotes on 31. March 2014."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The SRJD Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given is a filtered probability space $\\{\\Omega,\\mathcal{F},\\mathbb{F},P\\}$ representing uncertainty in the model economy with final date $T$ where $0<T<\\infty$. $\\Omega$ denotes the continuous state space, $\\mathcal{F}$ a $\\sigma-$algebra, $\\mathbb{F}$ a filtration &mdash; i.e. a family of non-decreasing $\\sigma-$algebras $\\mathbb{F}\\equiv\\{\\mathcal{F}_{t\\in [0,T]}\\}$ with $\\mathcal{F}_{0}\\equiv\\{\\emptyset,\\Omega\\}$ and $\\mathcal{F}_{T}\\equiv \\mathcal{F}$ &mdash; and $P$ the real or objective probability measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the SRJD model, which is an affine jump diffusion (see Duffie et al. (2000)), the risk-neutral dynamics of the VSTOXX volatility index are given by the following stochastic differential equation (SDE)\n",
    "\n",
    "\n",
    "$$\\begin{equation*}\n",
    "dv_{t}=\\kappa (\\theta -v_{t})dt+\\sigma \\sqrt{v_{t}}dZ_{t} + J_{t} v_{t} dN_{t}  -r_{J}dt\n",
    "\\end{equation*}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The meaning of the variables and parameters is\n",
    "\n",
    "* $v_{t}$ volatility index level at date $t$\n",
    "* $\\kappa$ speed of adjustment of $v_{t}$ to ...\n",
    "* ... $\\theta$, the long-term mean of the index\n",
    "* $\\sigma$ volatility coefficient of the index level\n",
    "* $Z_{t}$ standard Brownian motion\n",
    "* $J_{t}$ jump at date $t$ with distribution ...\n",
    "* ... $\\log(1+J_{t})\\approx \\mathbf{N}\\left(\\log(1+\\mu)-\\frac{\\delta^{2}}{2},\\delta^{2}\\right)$\n",
    "* $\\mathbf{N}$ cumulative distribution function of a standard normal random variable\n",
    "* $N_{t}$ Poisson process with intensity $\\lambda$\n",
    "* $r_{J}\\equiv \\lambda \\cdot \\left(e^{\\mu + \\delta^{2}/2}-1\\right)$ drift correction for jump"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stochastic process for $v$ is adapted to the filtration $\\mathbb{F}$. Moreover, $Z$ and $N$ are not correlated. The time $t$ value of a zero-coupon bond paying one unit of currency at $T,0\\leq t < T,$ is $B_{t}(T) = e^{-r(T-t)}$ with $r \\geq 0$ the constant risk-less short rate. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the Fundamental Theorem of Asset Pricing, the time $t$ value of an attainable, $\\mathcal{F}_{T}-$measurable contingent claim $V_{T}\\equiv h_{T}(X_{T})\\geq 0$ (satisfying suitable integrability conditions) is given by arbitrage as\n",
    "\n",
    "\\begin{equation*}\n",
    "V_{t}=\\mathbf{E}^{Q}_{t}\\left( B_{t}(T) V_{T} \\right)\n",
    "\\end{equation*}\n",
    "\n",
    "with $V_{0}=\\mathbf{E}^{Q}_{0}\\left( B_{0}(T) V_{T} \\right)$ as the important special case for valuation purposes. $Q$ is a $P-$equivalent martingale measure. The contingent claim could be a European call option maturing at $T$ with payoff $V_T = h_{T}(v_{T})\\equiv\\max[v_{T}-K,0]$. It could also be a European put with payoff $V_T = h_{T}(v_{T})\\equiv\\max[K-v_{T},0]$. In both cases, $K$ is the fixed strike price of the option. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate the financial model, i.e. to generate numerical values for $v_{t}$, it has to be discretized. To this end, divide the given time interval $[0,T]$ in equidistant sub-intervals $\\Delta t$ such that now $t \\in \\{0,\\Delta t,2 \\Delta t, ..., T\\}$, i.e. there are $M+1$ points in time with $M\\equiv T/\\Delta t$. With $s=t-\\Delta t$, a discretization of the continuous time market model is given by \n",
    "\n",
    "$$\\begin{eqnarray*}\n",
    "\\tilde{v}_{t} &=& \\tilde{v}_{s} + \\kappa(\\theta-\\tilde{v}_{s}^{+}) \\Delta t + \\sigma\\sqrt{\\tilde{v}_{s}^{+}} \\sqrt{\\Delta t}z^1_{t} \\nonumber \\\\\n",
    "&+& \\left(e^{\\mu_{J}+\\delta^{2}z^{2}_{t}}-1\\right) \\tilde{v}_{s}^+ y_{t} -r_J \\Delta t \\\\\n",
    "v_{t} &=& \\tilde{v}_{t}^{+} \n",
    "\\end{eqnarray*}$$\n",
    "\n",
    "for $t \\in \\{\\Delta t, ..., T\\}$ with $x^+\\equiv\\max[x,0]$ and the $z^{n}_{t}$ being standard normally distributed and $y_t$ Poisson distributed. $z^{1}_{t}, z^{2}_{t}$ and $y_t$ are uncorrelated. This discretization scheme is an Euler discretization and is generally called *full trunction* scheme. See Lord et al. (2008) for an analysis of this and other biased discretization schemes for the square-root diffusion process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Term Structure Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first in the  calibration of the SRJD model, is with regard to the futures term structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Futures Term Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is difficult for parsimonious short rate models like Cox et al. (1985) to account for different term structures of the interest rate. A possible solution is the introduction of time-dependent parameters which, however, enlarges the number of parameters significantly, sacrificing at the same time the convenience of a limited number of economic parameters. Another solution is a *deterministic shift approach* according to Brigo and Mercurio (2001) which preserves the basic structure of the model with all its advantages and which nevertheless allows to better account for different term structures of the short rate.\n",
    "\n",
    "In this section, we transfer the deterministic shift approach for short rate models of Brigo and Mercurio (2001) to the SRJD model. Since the square-root diffusion and the jumps are not correlated, we can first apply the approach to the diffusion part and use this enhanced component later on in conjunction with the jump component.\n",
    "\n",
    "Before we present the theory, a look at the VSTOXX futures data first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "futures_quotes = h5['futures_data']\n",
    "futures_quotes['DATE'] = pd.to_datetime(futures_quotes['DATE'])\n",
    "futures_quotes['MATURITY'] = pd.to_datetime(futures_quotes['MATURITY'])\n",
    "futures_quotes.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure presents the futures quotes for all eight maturities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ax = futures_quotes.plot(x='MATURITY', y='PRICE',\n",
    "                         figsize=(10, 6), legend=False)\n",
    "futures_quotes.plot(x='MATURITY', y='PRICE', style='ro', ax=ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">VSTOXX futures quotes on 31. March 2014."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider for now the square-root diffusion volatility model of Gruenbichler and Longstaff (1996) as presented in the previous chapter, which is formally the same as the short rate model of Cox et al. (1985) and which is a special case of the SRJD model\n",
    "\n",
    "\\begin{equation*}\n",
    "dv_t = \\kappa (\\theta -v_{t})dt + \\sigma \\sqrt{v_{t}}dZ_{t}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to calibrate this model to the observed volatility term structure, given by the above presented set of prices for futures on the VSTOXX index with different maturities. We have to minimize for all considered times $t$ and a parameter set $\\alpha = (\\kappa,\\theta,\\sigma,v_{0})$ simultaneously the single differences\n",
    "\n",
    "\\begin{equation*}\n",
    "\\Delta f(0,t) \\equiv f(0,t)-f^{GL96}(0,t;\\alpha)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $f(0,t)$ is the time $0$ market (instantaneous) forward volatility for time $t$ and the quantity $f^{GL96}(0,t;\\alpha)$ is the model (instantaneous) forward volatility for time $t$ given parameter set $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that there is a continuously differentiable volatility term structure function $F(0,t)$ available (i.e. there are infinitely many volatility futures prices). The forward volatility then is\n",
    "\n",
    "\\begin{equation*}\n",
    "f(0,t)=\\frac{\\partial F(0,t)}{\\partial t}\n",
    "\\end{equation*}\n",
    "\n",
    "On the other hand, the model implied forward volatility is given as (see Brigo and Mercurio (2001))\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "f^{GL96}(0,t;\\alpha) = \\frac{\\kappa \\theta\\left(e^{\\gamma t}-1\\right)}{2\\gamma+(\\kappa+\\gamma)\\left(e^{\\gamma t}-1\\right)} \\nonumber \\\\\n",
    "+v_{0} \\frac{4 \\gamma ^{2} e^{\\gamma t}}{\\left(2\\gamma+(\\kappa+\\gamma)\\left(e^{\\gamma t}-1\\right)\\right)^{2}} \n",
    "\\end{eqnarray*}\n",
    "\n",
    "with\n",
    "\n",
    "\\begin{equation*}\n",
    "\\gamma \\equiv \\sqrt{\\kappa^{2}+2\\sigma^{2}}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python script ``srjd_fwd_calibration.py`` contains the Python code to calibrate the forward volatilities to the VSTOXX futures prices (see the appendix for the complete script). The beginning of the script is about library imports, importing the data sets and making some selections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 8,28p scripts/srjd_fwd_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``srd_forwards()`` implements the forward formula from above for a given parameter set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('scripts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import srjd_fwd_calibration as srjdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdf.srd_forwards??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To operationalize the calibration, we use the mean squared error (MSE) as our yardstick\n",
    "\n",
    "\\begin{equation*}\n",
    "\\min_{\\alpha } \\frac{1}{N} \\sum_{n=1}^{N}\\left( f_{n} - f_{n}^{GL96}(\\alpha )\\right)^{2}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is to be minimized. Here, we assume that we have $N$ observations for the forward volatility. In Python this takes on the form of function ``srd_fwd_error()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdf.srd_fwd_error??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Executing the script yields optimal parameters for the Gruenbichler and Longstaff (1996) model given the VSTOXX futures prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%run scripts/srjd_fwd_calibration.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "opt.round(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These optimal values can be used to calculate the model forward volatilities using function ``srd_forwards()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from srjd_fwd_calibration import *\n",
    "srd_fwds = srd_forwards(opt)\n",
    "srd_fwds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numerical differences to the market futures prices are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "srd_fwds - forwards"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure compares the model futures prices (forward volatilities) with the VSTOXX futures market quotes. For longer maturities the fit is quite well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6));\n",
    "plt.plot(ttms, forwards, 'b', label='market prices');\n",
    "plt.plot(ttms, srd_fwds, 'ro', label='model prices');\n",
    "plt.legend(loc=0);\n",
    "plt.xlabel('time-to-maturity (year fractions)');\n",
    "plt.ylabel('forward volatilities');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">VSTOXX futures market quotes vs. model prices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we save the results from the term structure calibration for later use during the simulation of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "f = open('varphi', 'wb')  # open file on disk\n",
    "## write ttms object and differences (varphi values) as dictionary\n",
    "pickle.dump({'ttms': ttms, 'varphi': srd_fwds - forwards}, f)\n",
    "f.close()  # close file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shifted Volatility Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that we are given a continuously differentiable futures price function (i.e. through splines interpolation of discrete futures prices for different maturities). We consider now the deterministically shifted volatility process (see Brigo and Mercurio (2001))\n",
    "\n",
    "\\begin{equation*}\n",
    "\\hat{v}_t \\equiv v_t + \\varphi(t,\\alpha^*)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $\\varphi(t,\\alpha^*) \\equiv f(0,t)-f^{GL96}(0,t;\\alpha^*)$, the difference at time $t$ between the market implied forward volatility and the model implied forward volatility after calibration, that is for the optimal parameter set $\\alpha^*$. $\\varphi(t,\\alpha^*)$ corresponds to the differences (bars) in the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6));\n",
    "plt.bar(ttms, srd_fwds - forwards,\n",
    "            width=0.05, label='$\\\\varphi(t,\\\\alpha^*)$')\n",
    "plt.legend(loc=0)\n",
    "plt.xlabel('time-to-maturity (year fractions)')\n",
    "plt.ylabel('deterministic shift');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Deterministic shift values to account for VSTOXX futures term structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SRJD model discretization can now be adjusted as follows:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\tilde{v}_{t} &=& \\tilde{v}_{s} + \\kappa(\\theta-\\tilde{v}_{s}^{+}) \\Delta t + \\sigma\\sqrt{\\tilde{v}_{s}^{+}} \\sqrt{\\Delta t}z^1_{t} \\nonumber \\\\\n",
    "&+& \\left(e^{\\mu_{J}+\\delta^{2}z^{2}_{t}}-1\\right)  \\tilde{v}_{s}^+ y_{t} -r_J \\Delta t  \\label{eq:disc3} \\\\\n",
    "\\hat{v}_{t} &=& \\tilde{v}_{t}^{+} + \\varphi(t,\\alpha^*) \\label{eq:disc4}\n",
    "\\end{eqnarray}\n",
    "\n",
    "This is consistent since the diffusion and jump parts are not correlated and since the jump part is added in a way that the first moment of the stochastic volatility process does not change."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Option Valuation by Monte Carlo Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section implements Monte Carlo simulation procedures for the SRJD model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Valuation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In what follows, the model option values are computed by Monte Carlo simulation (MCS). Given the discrete version of the financial model, the value of a European call option on the volatility index is estimated by MCS as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://hilpisch.com/lvvd_algo.png\" align=\"center\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$V_0$ is the MCS estimator for the European call option value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Technical Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In view of using a numerical method like MCS for all valuation and calibration tasks, the parametrization and implementation of the MCS algorithm plays an important role. Some major features of our implementation are:\n",
    "\n",
    "* **discretization**: the algorithm uses the Euler discretization scheme which is an approximation only but which might bring performance benefits\n",
    "* **random numbers**: for every single option valuation the seed can be held constant such that every option is valued with the same set of (pseudo-) random numbers\n",
    "* **variance reduction**: both antithetic variates and moment matching (for the first two moments of the pseudo-random numbers) are used as generic variance reduction techniques\n",
    "* **deterministic shift**: the deterministic shift values $\\varphi$ are determined only once through a separate calibration and are held constant afterwards (even if model parameters change)\n",
    "\n",
    "$\\varphi$ only has to be deterministic and integrable on closed intervals (see Brigo and Mercurio (2001)), which is of course the case. For the approach to be valid, it is not important how we came up with the $\\varphi$ originally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python code for simulating the SRJD model is found in the script ``srjd_simulation.py`` (see the appendix for the complete script). The beginning of the script shows several imports, the definition of example parameters and also the cubic splines interpolation to be used for the estimation of the deterministic shift parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 8,39p scripts/srjd_simulation.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Python function ``random_number_gen()`` generates arrays of standard normally distributed pseudo-random numbers using both antithetic variates and moment matching as generic variance reduction techniques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import srjd_simulation as srjds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjds.random_number_gen??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The major function of this script is ``srjd_simulation()`` which implements the Monte Carlo simulation for the SRJD model based on an Euler discretization scheme. The scheme used here is usually called *full truncation* scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjds.srjd_simulation??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the function ``srjd_call_valuation()`` estimates the value of a European call option given the simulated volatility paths from ``srjd_simulation()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjds.srjd_call_valuation??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Executing the script yields a MCS estimator for the European call option with the parameters as assumed as in the script of about 1 currency unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%run scripts/srjd_simulation.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section now calibrates the SRJD model to market quotes for European call options on VSTOXX futures. It considers calibrations to a single maturity as well as to multiple maturities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Python Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calibration of the SRJD model is similar to the procedure for the Gruenbichler and Lonstaff (1996) square-root diffusion model as presented in the previous chapter. The major difference now is that we have to take into account more parameter values for the optimization. The Python code is contained in script ``srjd_model_calibration.py`` (see the appendix for the complete script). As usual a few imports and parameter definitions first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 10,21p scripts/srjd_model_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In what follows, we want to calibrate the model simultaneously to multiple maturities for the VSTOXX European call options. The valuation function ``srjd_valuation_function()`` therefore now calculates the differences between model and market values directly and returns an array with all differences (relative or absolute)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import srjd_model_calibration as srjdc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdc.srjd_valuation_function??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error function ``srjd_error_function()`` for the SRJD model has to be enhanced compared to the SRD case to account for the additional parameters of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdc.srjd_error_function??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same holds true for the calibration function ``srjd_model_calibration()`` itself. This function allows to select certain maturities for the the calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdc.srjd_model_calibration??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Short Maturity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The addition of a jump component shall allow a better fit to short term call option market quotes. Therefore, consider the following calibration to the shortest option maturity available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from srjd_model_calibration import *\n",
    "## read option data, allow for 30% moneyness tolerance\n",
    "option_data = read_select_quotes(tol=0.3)\n",
    "option_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "opt_1 = srjd_model_calibration(option_data, p0=None,\n",
    "                               rel=False, mats=['2014-4-18'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal parameter values are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "opt_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these optimal parameter values, add the model prices to the ``DataFrame`` object containing the option data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "values = []\n",
    "kappa, theta, sigma, lamb, mu, delta = opt_1\n",
    "for i, option in option_data.iterrows():\n",
    "    value = srjd_call_valuation(v0, kappa, theta, sigma,\n",
    "                                lamb, mu, delta,\n",
    "                                option['TTM'], r, option['STRIKE'],\n",
    "                                M=M, I=I, fixed_seed=True)\n",
    "    values.append(value)\n",
    "option_data['MODEL'] = values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the calibration results graphically. Indeed, the fit seems to be quite good, reflecting a MSAE of about 0.0015 only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "## selecting the data for the shortest maturity\n",
    "os = option_data[option_data.MATURITY == '2014-4-18']\n",
    "## selecting corresponding strike prices\n",
    "strikes = os.STRIKE.values\n",
    "## comparing the model prices with the market quotes\n",
    "fig, ax = plt.subplots(2, 1, sharex=True, figsize=(10, 6));\n",
    "ax[0].plot(strikes, os.PRICE.values, label='market quotes');\n",
    "ax[0].plot(strikes, os.MODEL.values, 'ro', label='model prices');\n",
    "ax[0].legend()\n",
    "ax[1].bar(strikes, os.MODEL.values - os.PRICE.values,\n",
    "                        width=0.3);\n",
    "ax[1].set_xlim(12.5, 23);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of SRJD model to European call options on the VSTOXX for shortest option maturity (April 2014)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two Maturities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us proceed with the simultaneous calibration of the model to the May and July maturities. The MSE is also pretty low in this case (i.e. below 0.01)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "## read option data, allow for 17.5% moneyness tolerance\n",
    "option_data = read_select_quotes(tol=0.175)\n",
    "# option_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "opt_2 = srjd_model_calibration(option_data, rel=False,\n",
    "                            mats=['2014-5-16', '2014-7-18'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal parameter values are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "opt_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In what follows, we use the Python function ``plot_calibration_results()`` to generate the plots for the different valuation runs. This function allows for different numbers of sub-plots, i.e. when the number of option maturities is changed. It is mainly a generalization of the plotting code used above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srjdc.plot_calibration_results??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the results of the  calibration graphically. It is obvious that the SRJD model is able to account for multiple maturities at the same time which is mainly due to the term structure component introduced by the deterministic shift approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(option_data, opt_2, ['2014-5-16', '2014-7-18'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of SRJD model to European call options on the VSTOXX for May and July 2014 maturities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Four Maturities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a next step, we consider four maturities &mdash; before we try to calibrate the model to all eight maturities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "mats = sorted(set(option_data['MATURITY']))\n",
    "mats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The four maturities for this particular calibration run are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "mats[::2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this calibration run, we use the optimal parameters from the previous calibration to two maturities. Obviously, the more options to calibrate the model to, the longer the procedure takes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "opt_4 = srjd_model_calibration(option_data, p0=opt_2,\n",
    "                               rel=False, mats=mats[::2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "opt_4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even calibrating the model to four maturities yields quite a good fit over these maturities as the following figure illustrates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(option_data, opt_4, mats[::2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of SRJD model to European call options on the VSTOXX for four maturities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All Maturities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let us attack the hardest calibration problem &mdash; the one involving all eight option maturities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "opt_8_MSAE = srjd_model_calibration(option_data,\n",
    "                                rel=False, mats=mats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "opt_8_MSAE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following shows that the fit is still reasonable for eight maturities and that many options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(option_data, opt_8_MSAE, mats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of SRJD model to European call options on the VSTOXX for all eight maturities (MSAE used)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check whether there is a (larger) difference when we calibrate the model using relative differences (i.e. the MSRE) as yardstick, consider the following calibration run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "opt_8_MSRE = srjd_model_calibration(option_data, p0=opt_8_MSAE,\n",
    "                                    rel=True, mats=mats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure presents the results. They are not too dissimular to the ones obtained using the MSAE as yardstick. The major difference is the weighting of the options in that now those options with lower market quotes (higher strikes) get more weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plot_calibration_results(option_data, opt_8_MSRE, mats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style=\"font-family: monospace;\">Calibration of SRJD model to European call options on the VSTOXX for all eight maturities (MSRE used)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter introduces a more sophisticated model, the so-called square-root jump diffusion model (SRJD), for the evolution of the VSTOXX volatility index over time. It enhances the Gruenbichler and Longstaff (1996) square-root diffusion model by two components: a log normally distributed *jump component* and a *determistic shift component*. While the first allows to better calibrate the model to short term option quotes, the latter makes it possible to take the volatility term structure &mdash; as embodied by the eight futures on the VSTOXX index &mdash; into account. All in all, the model yields good calibration results even in cases where all eight option maturities are accounted for."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Scripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### srjd_fwd_calibration.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srjd_fwd_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### srjd_simulation.py "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srjd_simulation.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### srjd_model_calibration.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat scripts/srjd_model_calibration.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://hilpisch.com/tpq_logo.png\" alt=\"The Python Quants\" width=\"35%\" align=\"right\" border=\"0\"><br>\n",
    "\n",
    "<a href=\"http://tpq.io\" target=\"_blank\">http://tpq.io</a> | <a href=\"http://twitter.com/dyjh\" target=\"_blank\">@dyjh</a> | <a href=\"mailto:team@tpq.io\">team@tpq.io</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
